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Abstract 

The objective of quantitative photoacoustic tomography (QPAT) is to reconstruct 
optical and thermodynamic properties of heterogeneous media from data of absorbed 
energy distribution inside the media. There have been extensive theoretical and com- 
putational studies on the inverse problem in QPAT, however, mostly in the diffusive 
regime. We present in this work some numerical reconstruction algorithms for multi- 
source QPAT in the radiative transport regime with energy data collected at either 
single or multiple wavelengths. We show that when the medium to be probed is non- 
scattering, explicit reconstruction schemes can be derived to reconstruct the absorption 
and the Griineisen coefficients. When data at multiple wavelengths are utilized, we 
can reconstruct simultaneously the absorption, scattering and Griineisen coefficients. 
We show by numerical simulations that the reconstructions are stable. 

Key words. Quantitative photoacoustic tomography (QPAT), sectional photoacoustic tomogra- 
phy, radiative transport equation, inverse transport problem, interior data. Born approximation, 
iterative reconstruction. 

1 Introduction 

In photoacoustic tomography (PAT) experiment, we send near infra-red (NIR) light into a 
biological tissue. The tissue absorbs part of the incoming light and heats up due to the 
absorbed energy. The heating then results in expansions of the tissue and the expansion 
generates compressive (acoustic) waves. We then measure the time-dependent acoustic signal 
that arrives on the surface of the tissue. From the knowledge of these acoustic measurements, 
we are interested in reconstructing the absorption and scattering properties of the tissue, as 
well as the thermodynamic Griineisen parameter which measures the photoacoustic efficiency 
of the tissue. We refer interested readers to [4, 12, 22, 28, 33, 59, 61, 67, 72, 95, 96, 100] for 
overviews of the field of photoacoustic imaging. 
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The propagation of NIR light in biological tissues is accurately modeled by the radiative 
transport equation which describes the distribution of photons in the phase space. To be 
precise, let Q ^ (d > 2) be the domain of interest with smooth boundary dQ and 
be the unit sphere in M^. We denote hy X = Q x S>^~^ the phase space and r± = {(x, v) : 
(x, v) ednx §^-1 s.t. ± n(x) • V > 0} its incoming and outgoing boundaries, n(x) being the 
unit outer normal vector at x G dQ. The radiative transport equation for photon density 
2i(x, v) can then be written as [J, 11, 77]: 

V • V2i(x, v) + cra(x)2i(x, v) = cr5(x)X(2i)(x, v) in X 
t^(x, v) = 5'(x, v) on r_. 

Here the positive functions cra(x) and as{'x.) are the absorption and the scattering coefficients, 
respectively. The function ^'(x, v) is the incoming illumination source. The scattering oper- 
ator K is defined as 

X(^)(x,v)= / /C(v,vXx,vOrfV-t^(x,v) 

where dv is the normalized measure on S'^-^ in the sense that J^^.i rfv = 1, and the kernel 
/C(v,v') describes the way that photons traveling in direction getting scattered into 
direction v, and satisfies the normalization condition /C(v, v')rfv' = 1, V v G S>^~^. 
In practical applications in biomedical optics, JC is often taken to be the Henyey-Greenstein 
phase function [9, 54, 97] which depends only on the product v • v'; see equation (5.5) in 
Section 5. 

The photon energy that is absorbed at location x G per unit volume, £'(x), is the 
product of the absorption coefficient and the ffuence distribution: 

E{^)= [ a,(x)^(x,v)rfv. (1.2) 

The heating due to this absorbed energy generates an initial pressure field, denoted by i7, 
in the tissue that depends on the thermodynamic property of tissue and is proportional to 
E: 

//(x) = T(x)£;(x) = T(x) / a,(x)^(x, v)rfv (1.3) 

where the positive function T(x) is the nondimensional Griineisen coefficient which in the 
current formulation measures the photoacoustic efficiency of the tissue. To simplify the 
presentation, we use the short notation (/)v to denote the integral of / over the v variable 
in the rest of this work. 

The initial pressure field H then propagates according to the acoustic wave equation, 
with the wave speed c(x), 

Ap = 0, in M+ xM^ 

c2(x) dt^ (14) 

j9(0, x) = //(x), ^(0, x) = in M^. 

The time-dependent pressure signal p{t^ x) is then measured on the surface of the tissue for 
long enough time, say t^o^ and the objective is to reconstruct the coefficients a^, cFs and the 
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Griineisen coefficient T from this measurement. Note that the reason that we can write the 
transport equation in stationary case while using the wave equation in time-dependent case, 
is that the two phenomena occur on two significantly different time scales [ ] . 

The reconstruction problem in photoacoustic tomography can be split into two steps. 
In the first step, we need to reconstruct the initial pressure field i/(x) from measured 
acoustic signal on the boundary, This is a well known inverse problem 

for the acoustic wave equation that has been thoroughly studied in the past a few years 
under various scenarios; see for instance [ , 24, 29, 40, 42, 4i, 4b>, bU, ^i, 52, o5, bU, 
64, 65, 62, 63, 68, 69, 72, 82, 99] for analytical reconstruction formulas with constant wave 
speed, [55, 56, 73, 85, 88, 89, 98] for reconstruction under variable wave speed, and [7, 8, 
36, 58, 91, 92, 3, 23, 32, 36, 71, 70, 73, 87, 86, 90, 101] for reconstructions under even more 
complicated environments. 

This work is concerned with the second step of photoacoustic tomography, call quanti- 
tative photoacoustic tomography (QPAT). The objective is to reconstruct the absorption 
and the diffusion coefficients, a a and a^, in the transport equation (1.1) and the Griineisen 
coefficient T from the result of the first step, the data H in (1.3). This step has recently 
attracted significant attention from both mathematical [5, 6, 13, 21, 14, 15, 16, 18, 19], 
computational [20, 30, 31, 34, 48, 44, bb, ^0, 81, 83, 104, 106] and modeling and experimen- 
tal Y ^-i 75] perspectives. Most of existing work on this step, however, is done in the diffusive 
regime [14, 16, 18], i.e., based on the diffusion approximation to the transport equation (1.1). 
Transport-based QPAT is only studied in [13, 34, 102] where the Griineisen coefficient T 
has been assumed to be a constant and known. 

To setup the problem appropriately, for the rest of the paper, we assume that all the 
coefficients that we are interested in are positive and bounded functions. More precisely, we 
assume: 

(Al) the coefficients < Cq < T(x), cra(x), ^^(x) G L^{Q) n L'^{Q) for some Cq > 0; 
(A2) the scattering kernel /C(v, V) G Li(§^-i x S^-^); 

(A3) the illumination, modeled by the boundary condition is in L^(r_, rf^), with measure 
(i^ = |v • n(x)|rfm(x)rfv, rfm(x) being the usual Lebesgue measure on d^. 

With these assumptions, it is well-known that the radiative transport problem (1.1) is well- 
posed and thus admits a unique solution u G L^{X) [ ]. This means that the data H in (1.3) 
is a well-defined function in L^(f]). In practical applications, the objects to be imaged are 
often embedded into phantoms of similar optical and acoustic properties to get regularly 
shaped imaging domains. We thus assume will assume that the domain Vt is convex. This 
assumption will simplify some of the presentation but is not essential for the results obtained. 

We conclude this section with the two remarks. First, if both the absorption coefficient 
aa(x) and the scattering coefficient ^^(x) are known and only the Griineisen coefficient has 
to be reconstructed, we can simply solve the transport equation (1.1) and compute the 
energy £'(x). Then T(x) is reconstructed as T = ^. Thus we need only one interior data 
set and one transport solver to solve the inverse problem. This is a trivial case. We will not 
discuss this case in the rest of the paper. 
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Second, in practical application of PAT in biomedical imaging, the absorption and the 
scattering coefficients and (jg are often isotropic, i.e. independent of the angular variable 
V. We thus restrict ourselves to the case of isotropic coefficients in this work. Mathematically 
this assumption is essential, as we can see from the following result that it is not possible to 
reconstruct uniquely anisotropic coefficients. 

Proposition 1.1. Let (T, a^, (Js) ci'^d (T, a^, a^) he two sets of coefficients, and H and H the 

corresponding data sets. Let z(x) G C^{^) he an arhitrary positive function with houndary 
value z\q^ = 1. Then 

5"a = (<^a — V • V In z)z and Ta^ = Ta^ (1.5) 

implies H = H. 

Proof Let u be the solution of the transport equation for coefficients (a^, CFg) with boundary 
value g. It is straightforward to verify that uz is the solution of the transport equation for 
coefficients (cr^ — v • V In z^ag) with the same boundary value g (because z\q^ = 1). Now, 

clearly T / (a^ — v • V In z)uz dv = T GaU dv if (1.5) holds. □ 

The rest of the paper is structured as follows. We first study in Section 2 the inverse 
transport problem for non-scattering media for applications in quantitative sectional QPAT. 
We present some analytical reconstruction strategies in this setting. We then study in Sec- 
tion 3 the same inverse problem for scattering media. We linearize the nonlinear inverse 
problem using the tool of Born approximation and present some numerical procedures to 
solve the linear and nonlinear inverse problems. In Section 4 we consider the QPAT prob- 
lem in the case when illuminations with multiple wavelengths are available. To validate the 
reconstruction strategies, we provide in Section 5 several numerical simulations with syn- 
thetic data generated under different scenarios. We conclude the paper with some additional 
remarks in Section 6. 



2 QPAT of non-scattering media 

We start by considering the QPAT problem for non-scattering media. In this case, photons 
propagate in the media along straight trajectories, without changing directions of propaga- 
tion. The scattering mechanism in the transport equation (1.1) is thus dropped by setting 
the scattering coefficient = 0. We have the following transport model for light propaga- 
tion: 

V • V2i(x, v) + cra(x)t^(x, v) = in X 

2i(x,v) = 5'(x, v) on r_. 

The fact that photons travel in straight lines allows us to illuminate only part of the domain, 
for instance a plane cut through a three-dimensional medium. This is the fundamental 
principle of sectional photoacoustic tomography [107, 39]. 

To simplify the presentation, for any point x G and direction v G let us define 

r±(x, v) = inf{s G ]R+|x± sv ^ f]}. 
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and r(x, v) = r+(x,v) + r_(x, v). It is easy to see that r+(x, v) (resp. r_(x, v)) is the 
distance it takes a particle at x travehng in direction v (resp. — v) to reach the boundary 
of the domain. In the same spirit, for any point (x,v) on the incoming boundary r_, we 
define 

r+(x, v) = sup{5 G ]R+|x ± 5V G f^}, 

and set r_(x, v) = 0. Thus r+(x, v) is the distance for a photon coming into the domain at 
X in direction v to exit the domain. 

It is straightforward to show the foUowing representation of the interior data. 

Lemma 2.1. The interior data H{x.) generated with source g and coefficients (T,^^) can 
be written as 

H{^) = T(x)a«(x) / ^(x - r_(x, v)v, v)e-^o^"^^'^^^^("-^-("'")^+^")^^rfv. (2.2) 

Proof. The transport equation (2.1) can be integrated along direction v as an ODE to obtain 
the solution: 

2i(x,v) = ^(x-r_(x,v)v,v)e-^o^"^^'^^^^("-^-("'")"+^")^^ (2.3) 
The result then follows from using this solution in data H introduced in (1.3). □ 

In fact, this simple representation of the transport solution allows us to obtain explicit 
procedure to reconstruct the two unknown coefficients T and if we have the luxury of 
using the right illumination source ^'(x, v). 




Figure 1: Two illumination schemes that allow the simultaneous reconstruction of T and 
aa using two data sets in a circular domain. Left: two collimated sources supported on S_ 
(solid part of dQ) and S+ (dashed part of dQ) respectively; Right: two point sources located 
at x' and x^' respectively. 



2.1 Reconstruction with collimated sources 

The first case where we can derive an analytical reconstruction formula is when collimated 
sources, i.e., sources focused on a specific direction, are used. Let be the direction that the 
source points to, and define S±(v') = {x G dQ\ ±n(x) • > 0}, S-(v') (resp. S+(v')) being 
the part of the boundary where lines in direction enters (resp. leaves) the domain 
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As mentioned in the Introduction, we assume that the domain Q is convex. In this case, 
each point x G S_(v') is uniquely mapped to a point y G S+(v') as: y = x + r+(x, v')v' 
and vise versa: x = y + r+(y, — v')(— v'). 

The colhmated iUumination source that we choose has the form 

^(x,v)=0(x)5(v-v'), xeE_(v'). (2.4) 
This leads to the following expression for the interior data: 



Let x' G dQ be the track back of x to the boundary in — v direction, i.e., x' = x — r_(x, v')v', 
then H can be written as 

i/(x' + r_(x,v')v') _ ,^(,. + ,_(,,V)v')e-/o-*''■^''-(-'+-')^^ (2.5) 



0(x')T(x' + T_(x,v')v') 



(i) Reconstruction of aa. If T is known, then we can take the log of the above formula 
and differentiate with respect to t_ to obtain the following equation for aa- 

v' ■ V(7„(x) - a^x) - f V ■ V In ^) (t„(x) = 0, m Q 

where g in the term • Vln— is evaluated at x — r_(x, v')v^ This is an initial value 

problem for a first order ordinary differential equation (because the direction is fixed). 
It can be solved uniquely [ ] to reconstruct the absorption coefficient cr^ along the lines in 
direction v^ 

There is an equivalent reconstruction procedure for (7^. When a collimated source is 
used, photons travel only along the direction that the source is pointing to, say v^ Thus, the 
solution u to the transport equation is only non-zero in this direction. In other words, {u)^ = 
/gd-i '^(x, v)rfv = 2i(x, v'), and hence the data can be expressed as i7 = Taa(x)2i(x, v'). We 
can thus replace the absorption term aa(x)t^(x, v) in the transport equation with the term 
H/T to obtain the following transport equation for the photon density u: 

V • V2i(x, v) + Y = in X 7) 

2i(x, v) = g(x)5(v — v') on r_. 

This transport equation can be solved uniquely for u [ . .]: 

«(x' + tv, v) = I S(^') + ^ = (2.8) 

( 0, x' e E_ and v ^ v'. 

Once u is obtained, we can reconstruct the absorption aa = H/ (Tu), provided that we select 
the boundary condition g such that the denominator does not vanish. 
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(ii) Reconstruction of (T, a^). In fact, we can use two data sets generated with coUimated 
sources to determine the Griineisen coefficient and the absorption coefficients simultaneously. 
The configuration of the two sources is illustrated in Fig. 1 (left). Let 5'i(x, v) = g{x.)5{-v—V) 
be the first coUimated source supported on S_, and Hi be the corresponding interior data 
which takes the form of H in (2.5). Then the second source we use is supported on S+ 
with identical intensity distribution but pointing in the opposite direction, i.e., — v^ More 
precisely, 

^2(x,v)=0(x-r+(x,v')v')<5(v + v'), xeE+(v'). (2.9) 

If we use this source in the expression for the data H in (2.2), we obtain the following 
expression for the interior data H2: 

.'^mtr'r''T?. = + r_(x,v')v>- (2.10) 

0(x')T(x' + r_(x, v')v') 

We now take the logarithm of the ratio of (2.5) and (2.10) and use the relation r(x, v) = 
r+(x, v) + r_(x, v) to obtain 

TT /«r+(x^vO /'r_(x,vO 

ln--^(x' + r_(x,v>0 = - / aa{^' + sV)ds + 2 / aa{^' + sV)ds. (2.11) 

Hi Jo Jo 

Differentiation of this result with respect to r_(x, v') (equivalent to the directional differen- 
tiation • Vx) will allow us to reconstruct the quantity 2aa(x' + r_(x, v')v') = 2cra(x) a.e. . 
The coefficient T can then be reconstructed by solving the transport equation with the recon- 
structed aa and gi (resp. 5^2)5 and computing T = Hi/{aa{ui)^) (resp. T = i?2/(<^a('^2)v))- 

2.2 Reconstruction with point sources 

The second case where we can derive an analytical reconstruction method is when the 
illumination source is a point source in the spatial variable: 

p(x,v) = 0(v)5(x-x'). (2.12) 



X — x' 



This is the most commonly used type of source in optical imaging, such as optical tomogra- 
phy [9]. For any point x G f^, let us define 
to X. Then the interior data H simplifies to 



phy [9]. For any point x G f^, let us define = r as the unit vector pointing from x' 

X — x' 



i/(x) = 0(vO(T(x)a,(x)e-^o^"^^'^'^"^("-"-("'"')"'+^"')^^). 
Following the presentation in the previous section, we can rewrite this as 

"}r^":^Y^^'l, = + r.(..V)V)e-f;-"-''"^^'*'->". (2.13) 

0(vOT(x' + r_(x, v^v^ 

This is the same type of formula as (2.5) except that the unit direction vectors are different 
at different points x G f^. 
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(i) Reconstruction of a^. If T is known, we can differentiate (2.13) as before to obtain an 
equation for the absorption coefficient a^. For each G Sl~'^(x') = {v G S^~-^|n(x') • v < 0}, 

define Q^f = < x G 



|x-x^| 



>, then 



V . Va,(x) - a^(x) - • Vln Y^^) ^«(^) = 0' (2.14) 

Cra(x) = Cra(x'), at x'. 

For a fixed v', we can solve this first order ODE to find cra(x) along the line segment 
if we know the value of cra(x) at x^ Due to the singularity of the source at x', we can not 
reconstruct this value from data H as before. 

Point sources emit photons that travel in straight lines away from the source location. 
Thus the transport solution at each spatial position x is only non-zero in direction v^ Hence 
H = T(x)aa(x)i^(x, v'). The transport equation can again be simplified to: 

v-V2i(x,v) + Y = in X ^2.15) 

2i(x, v) = 0(v)5(x — x') on r_. 

This transport equation can be conveniently solved in the polar coordinates with the origin 
at x^ The absorption coefficient can then be reconstructed as = H/{Tu) away from x^ 

(ii) Reconstruction of (T, (7^). If both coefficients are unknown, we can use two data sets 
generated with point sources to reconstruct them. The setup is depicted in Fig. 1 (right). 
Let 5'i(x, v) = 0(v)(^(x — x') and 5^2 (x, v) = 0(v)(^(x — x^') be the two point sources used to 
produce the interior data. We denote by rf(x',x'') = |x' — x^''! the distance between the two 
x'' - x' 

points and by v = -— the unit vector pointing from x' to x'^ For any point x G fi. 



x-x'^ 



we define v^' = as the unit vector pointing from x^' to x. It is then straightforward 

|x — x^'l 

to verify that 

r_(x, vOV . V + r_(x, V^v'' • (-v) = rf(x^x'0• (2-16) 

The first source gi produces data Hi that is given in (2.13), while the second source g2 
produces the data H2 that is given by 

We can rewrite H2 into 

i/,(x" + r_(x,v")v") ^^(,. + ^_(,,V0v")e-/o-''''^"'-('^"+-")''^ (2.17) 



0(v")T(x" + r_(x,v")v") 
Taking the ratio of (2.17) and (2.13), we obtain, after taking logarithm on both sides 

^ , i/2(x" + r_(x, v")v")0(v')T(x' + r_(x, v')v') 
in 



Hi{-K' + r_(x, v')v')0(v")T(x" + t_(x, v")v") , 

r_(x,v') /-r-Cxy) 



rr-(:x.,v ) /.T_^x,v■■j 

/ aaix.' + sV)ds - / aaix." + sv")ds. (2.18) 

Jo Jo 
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We can now differentiate this result with respect to r_(x, v') (equivalent to the directional 
differentiation • Vx) and use relation (2.16) to obtain the quantity (1 — ^ 1 cra(x). 

Note that for any point x G located on the line determined by and x^', = — 1 
(because x locates between x' and x^'). The result in this case reduces to that in the case of 
two collimated sources: (1 — ^^^)aa{^) = 2cra(x). 



2.3 Reconstruction with cone-limited sources 

The reconstruction procedures in the previous sections work only when the illumination 
function g takes the required special forms. For more general sources, we do not have 
similar explicit reconstruction procedures. In fact, it is not clear whether or not data from 
an arbitrary source would uniquely determine the absorption coefficient. We consider here 
a reconstruction method for a source that is slightly more general than the sources used in 
Sections 2.1 and 2.2 but still possess certain causality properties of those sources. 

Let be a selected direction pointing inside the domain Vt. We intend to construct a 
source such that the transport equation (2.1) with the source is casual along direction v°. 
Such causality would allow us to derive a direct layer peeling method that solves the inverse 
problem in one pass in a non-iterative manner. 



(0, L) 



(U L) 




n 



<5<xxxxX' 



Figure 2: Domain J7 and its extension to a hyper-cube Vt. The physical sources located on 
a line marked with x are separated from by a hyper-plane IT. The boundary condition 
is non-zero only on the bottom side of Vt (thick black part) that is orthogonal to v^. The 
half-aperture 9^ of the cone V^(v^, 9^) is given by 9^ = 7r/2 — minj^i, ^2}- 
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We require the source ^'(x, v) satisfy the foUowing condition: 

^(x,v) = for any v ^ ^(v^ 0°) = {v G S^-^ v • v° > cos ^^j, (2.19) 

where the cone V{v^^9^) is determined by the selected direction G E>^~^ and the half- 
aperture < 9^ < 7t/2. a practically important case that satisfies the assumption (2.19) 
is when the physical source can be separated from the object occupying the domain Q 
by a hyper-plane. This limits the angle at which Q is visible from the source, thus the 
rays entering Q are limited to some cone of directions V{v^^9^) with orthogonal to the 
separating hyper-plane. Then the resulting effective boundary conditions satisfy (2.19). In 
fact, in such case we can extend to a rf-dimensional hyper-cube Cl such that ^'(x, v) is 
non-zero only on one side of the cube which is orthogonal to v^, as shown in Fig. 2. 

To simplify the presentation of the method we consider here the case = 2, although 
the algorithm remains essentially the same for = 3. We assume that Q is the square [0, L]^ 
with cra(x) = for X G We also assume that the rays enter the domain from the 

bottom side {{x^y) \ x G (0,L),y = 0}, so v° = (0, 1). The source ^'(x, v) is only non-zero 
on the bottom side. Thus the equation (2.1) is causal in the y direction; see Fig. 2 for the 
details of the setup. We use the parametrization v = (cos ^, sin ^) so that v G V{v^^9^) 
becomes 9 e [7r/2 - 9^ ,7r/2 + 9^]. This implies that we can divide (2.1) by sin^ (> 0) so 
that the free transport equation together with its boundary condition can be written as 

Uy{x,y,9) = -cote u^{x,y, 9) - csc9 aa{x,y) u{x,y, 9), 

u{x,o,9) = g{x,9), xe{o,L), ee[^-e',^ + e^] 

u{0,y,9) = 0, ye{0,L), ^ e [f - ^M) ^^'^^^ 
u{L,y,9) = 0, ye{0,L), ^G(f,f + ^'^ 



where we observe that no boundary condition is needed on the top side of the domain 
{(x, ^) I X G (0, L)^y = L}. The interior data i7(x, y) now takes the form 

H{x,y) = T{x,y)aa{x,y) J u{x,y,9)d9 = T{x,y)aa{x,y){u{x,y,9))eo. 

7r/2-6>o 

The role of the temporal variable in (2.20) is played by y (the system is causal in y) so we 
can apply a first order time stepping procedure to obtain a semi-discrete inversion scheme. 
Let us discretize (2.20) in ^ on a grid with nodes 7/^,^ = 0,..., A^^, where yo = 0, yNy = L 
and the grid steps are hk = yk — Vk-i^ = 1, . . . , A^^. For the semi-discrete quantities we 
use notation 

u^^\x,9)^u{x,yj,9), ai^\x)^aa{x,yj), j = 0,...,Ny. (2.21) 

To reconstruct with T known, we apply a forward Euler time stepping scheme to (2.20) 
to obtain the following explicit reconstruction procedure. 

1. Initialize u^^\x^ 9) = g{x^ 0), and 
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2. For j = 1, . . . ,Ny compute 

9) = {I- hj cot 9 8^- hj esc 9 ai^-^\x)) u^^-^\x, 9) (2.22) 
for X e (0, L) and 6 e [k/2 - 7r/2 + ^°]. At the boundary x e {0, L} set 



M(i)(0,^) = 0, ee [f -^°,f) 

M«)(L,^) = 0, ^G(f,f- 



(2.23) 



and compute the reconstruction 



The above method is quite stable in practice as demonstrated in the numerical examples 
in Section 5.1. However, the discretization of dx in (2.22) should be compatible with the 
grid refinement in y to retain stability of the forward Euler. If a coarser discretization in y 
compared to that in x is desired, then (2.22) can be replaced by a semi-implicit scheme 

(X + hj cot e dx) u^^\x, 9) = {I- hj 

which requires solving a first order differential (diflFerence if discretized) equation in x at 
every step j with boundary conditions (2.23). Generalizing this method to reconstruct 
simultaneously two coefficients T and remains a topic of future study. 

We conclude Section 2 by summarizing the reconstruction results for non-scattering 
media that we have introduced into the following uniqueness and stability theorem. 

Theorem 2.2. Let (T, a^) and (T, da) be two sets of coefficients in (1.3) and (2.1) satisfying 
the assumptions in (Al)^ and 5'i(x, v) and 5'2(x, v) be two source functions of the form (2.4) 
or (2.12). Let (i/i,i/2) ctnd {Hi^H2) be the corresponding interior data. Then (a) Hi = Hi 
implies ^^(x) = 5-a(x) ifT = T; and (b) [Hi,H2) = (^1,^2) implies (T,^^) = (T,5-a). 
Moreover, the following stability result holds: 

||T(x)a«(x)e-^o^"^^'^ - t(x)a«(x)e-^o^"^^'^ ^ ^^("-^-("'"')"'+^"')^^|Uoo(^) 

<Ci\\Hi-Hi\\loo^^), / = 1,2 (2.25) 

Ci being a constant that depends on gi but is independent of the data Hi and Hi . 



3 QPAT of scattering media 

We now consider the QPAT problem for scattering media. When the scattering effect is 
very weak, the results obtained in the previous section could be used to obtain a good 
approximation of the reconstruction. We thus assume that the scattering effects are suf- 
ficiently strong so that neglecting them would deteriorate significantly the quality of the 
reconstructions. On the other hand, we assume that the scattering effect is not strong 
enough for us to model the light propagation process with the diffusion model for which 
explicit reconstruction strategies have been proposed [14, 16, 18]. 
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3.1 Stability and uniqueness results 

Unlike the non-scattering regime where we have unique and stable reconstruction using 
only a small number of interior data sets, in scattering media, we have only stability and 
uniqueness results from information of the full albedo operator: 

It is the clear that A depends on all three coefficients T, and ag. We denote by 
||-4|U(Li(r_,dO;Li(f7)) the norm of ^ as a linear operator from L-^(r_,rf^) to L^{Q). 

We then have the following stability estimate following the analysis in [ , Theorem 2.3]. 

Theorem 3.1. Let (T, 0-^,^5) and (T,(7a,5"s) be two set of coefficients satisfying the reg- 
ularity assumptions in (Al)^ and A, A the corresponding albedo operators. Denote by 
a = aa + cFg and a = da + ^s- Then the following result holds 

^0 

< \\A- A\\c{L^{T_^d^).L^{p^)) (3.2) 

for any point (x^, v') G r_. 

The result is a slight modification of the result in [1 where the singular part of the 
Schwartz kernel of the albedo operator is analyzed in detail. The only change that we have 
to make is to include the Griineisen coefficient T in the final step of the analysis. We thus 
would not repeat the lengthy analysis here but refer to [13]. 

It turns out that the stability result (3.2) leads to the unique reconstruction of the 
coefficients in some simplified situations. More precisely, if one of the three coefficients is 
known, we can reconstruct the other two coefficients uniquely as summarized in the following 
corollary. 

Corollary 3.2. Let (T,^^,^^) and {T^da^cfs) be two set of coefficients satisfying the regu- 
larity assumptions in (Al). Then the following statements hold: (a) If T = T , then A = A 
implies (^^,^5) = (^^,5-5); (b) If as = then A = A implies (T,cr«) = {f,aa); (c) If 
= then A = A implies (T, as) = (T, 5"s). 

Proof Theorem 3.2 implies that we can reconstruct uniquely 

^^(x^v^^) = T(x' + ^vOa,(x' + ^vOe-^o"(-^+^vOrf^ (x^vO G r_. (3.3) 

We can reconstruct the same quantity for the a different point on (x' + r+(x', v')v', — v') on 
r_ with t replaced by r+(x', v') — t: A^(x' + r+(x', v')v', — v', r+(x', v') — t). Taking the log 
of the ratio of the two quant ites, we obtain 

M{^' + r+(x^ V)V, -V, r+(x^ V) - t) ^ 

M(x^v^^) 

/>r+(x^vO-^ rt 

- / a(x' + r+(x^ V) - sV)ds + / a(x' + sV)ds. (3.4) 

Jo Jo 
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Taking the derivative of the above quantity gives us 2cr(x' + tV) for a.e. (x^, v') G r_ and 
t > 0. This in turn gives us a(x) for a.e. x G f^. Once cr(x) is uniquely determined, we 
reconstruct Ta^ = ji uniquely from A^(x', V ^t). (a) If T is known already, then a a is known 
from /i, and thus a — a a would give us a^; (b) If cr^ is known, a — would give us 0-^. T 
is then known from //; (c) If a a is known, then ji would give us T and a would give us a^. 
This completes the proof. □ 

Remark 3.3. Case (a) was presented in [ ]. The construction of A^(x', v',t) and A^(x' + 
r+(x', v')v', — v', r+(x', v') — t) is in the same spirit as the construction of (2.5) and (2.10) 
(or (2.13) and (2.17)) using collimated (or point) sources in the previous section. 

The above stability and uniqueness results, however, do not guarantee the unique recon- 
struction of all three coefficients simultaneously. In fact, theory based on the diffusion model 
in [14] states that one can not simultaneously reconstruct all three coefficient uniquely unless 
additional information is available [ o]. Moreover, unlike the situation in the non-scattering 
regime, no explicit reconstruction method can be derived in the scattering media. We thus 
have to rely mainly on other computational methods for the reconstructions. We now de- 
rive a numerical reconstruction method based on the Born approximation of the nonlinear 
inverse problem. 



3.2 Linearized reconstruction with Born approximation 



We linearize around some known^ not necessarily constant, background optical properties 
aao(x) and aso{'x.). To be more precise, we assume 



cr«(x) = a«o(x) + 5-«(x). 



where the perturbations are small in the sense that 



a,(x) = a,o(x) + 5-5(x) 
5-a(x) 



(3.5) 



<C 1 and 



5-.(x) 



1. Note again that T(x) is not assumed in perturbative form since the inverse problem of 
reconstructing T is linear (because the data H is linearly proportional to the unknown T). 

The solution of the radiative transport problem with coefficients and can then be 
written as 

t^(x, v) = [/ (x, v) + ii(x, v), (3.6) 

where [/(x,v) is the solution of the transport equation (1.1) with the known background 
coefficients aao and a^o, and 'u(x, v) is the perturbation in the solution caused by the per- 
turbations 5-a(x) and 5-s(x) in the coefficient. The equation satisfied by the perturbation 
ii(x, v), to the first order, is 



V • V£t(x, v) + crao{-x.)u 
'u(x, v) 



aso{^)K{u) - da{^)U + ds{^)K{U) in X 
on r_ 



(3.7) 



It is can be shown, using the fact that u is Frechet differentiable with respective to and 
Gs [37, 38, 79], that the terms omitted are indeed high order terms. 



13 



We now introduce the (adjoint) Green's function G(x., v; y) for the transport problem 
with background optical properties, the solution of the following adjoint transport equation 



-V • VG(x,v;y) + crao(x)G'(x,v;y) = cr^o(x)/^(G)(x, v; y) - (5(x - y) inX 

G'(x,v;y) = on r+, 

(3.8) 

where S is the usual Dirac delta function. Note that since the transport operator is not 
self- adjoint, the boundary condition is now put on r+. Rigorously speaking, this equation 
only holds in the weak sense. We need to multiply it with a test function that is C+. That 
requires that the coefficients and a^o are at least continuous in Q. 

We can now multiply (3.7) by TG and integrate over X, and multiply (3.8) by Tu and 
integrate over X, and subtract the results to show that 

T(x)({i(x,v))v= / T(y)a„(y)([/(y,v)G(y,v;x)My 
Jn 

T(y)a,(y)(i^ (t/)(y, v)G(y, v; x)).dy. (3.9) 



The perturbations in (3.5) and (3.6) imply that the interior data H is now given, to the 
first order, by 

H = Ta,o(t/)v + Ta,([/)v + Ta,o(t^)v. (3.10) 
Combining (3.10) and (3.9), we can show that 

^ (x)=J(T) + £«(Ta,) + £^(Ta,), (3.11) 



where I is the identity operator and the operators and are defined, respectively, as 

^'^(T^a) = ^+ /T(y)a.(y)i^(y;x)dy (3.12) 

C%Ta,) = -^T(y)a.(y)^^^|^p^(y;x)rfy. (3.13) 

Here we have normalized the data by (Jao{U)v This can be done when ([/)v(x) 7^ which 
can be guaranteed by selecting appropriate illumination source 5'(x,v). The normalization 
makes the sizes of the kernels of the integral operators jC^ and jC^ on the same order at all 
locations. This is important due to the fact that the strength of optical signals usually varies 
over several orders of magnitude across the domain of interest. The normalization results in 
well-balanced matrix elements when the integral equation (3.11) is discretized for numerical 
solution. 

Equation (3.11) is a linear integral equation for the three variables T, Ta^ and Ta^. The 
kernels for the operator and that for {UG)^/ {U)^ and {K{U)G)^ / {U)^ are known 
since they only involve the solutions of the forward and adjoint transport equations with 
background optical properties 0-^0 and a^o. It remains to solve (3.11) to reconstruct the 
unknowns (T, Ta^ and Ta^), and then the real coefficients (T, da and dg). 



14 



In practice, we do not have the fuU albedo data, but only a finite number of data sets 
generated from different sources. Let us assume that there are A^^ data sets collected from 



A^^ illuminations. The data sets are denoted by {Hi}fj^ and the corresponding illuminations 



are denoted by {gi} 
written as 



No 



=1- 



The system of linear integral equation with A^^ data sets can be 








CI 














( T 


I 




















u 


ra 


rs 


) 





1 



Hi/{Ui)^ = hi 



(3.14) 



where and Cf are the same operators defined in (3.12) and (3.13) respectively, with U 
replaced with Ui. 

The integral formulation can be discretized to get a linear system of equations. Let us 
assume that we have a numerical procedure, say a quadrature rule, to discretize the integral 



equation (3.9), and assume that we discretize 5"a(x) on a mesh of A^^ nodes, {yk}k=i- Then 
C will be matrix consists of A^^ x 3 blocks, each having size A"^ x A"^. The Ik elements of 
the matrices are given by 



Ik 



Ik 



+ 6 



(^i(yfc, v)G'(y/c, v;x;)), 

(?7i(xz,v))v 



Ilk 



-6 



(K(^,)(yfc,v)G'(yfc,v;xO), 

([/i(xz, v))v 



(3.15) 

where ^ k < Nq) is the weight of the quadrature on the k-th element. For simplicity 

of notation, we will not differentiate from now on integral operators and their discrete 
equivalences, the matrices. 

The system (3.14) can be over- or under-determined and so is often solved in regularized 
least-square sense. More specifically, we solve the problem by solve the following minimiza- 
tion problem 

T \ / ^1 \ / T 



mm 

,Tcra,Tcrs 





|2 

1/2 



(3.16) 



where the first term is the data fidelity term and the second term is a Tikhonov regular- 
ization term with the strength of regularization given by p. The operator V is the discrete 
differentiation. The regularization term is needed due to the presence of noise in practice 
even though the inverse problem here is very stable compared to similar inverse transport 
problems in diffuse optical tomography [9, 11, 77]. 

Remark 3.4. The least squares formulation and Tikhonov regularization strategy (3.16) 
are adopted here mainly for their computational efficiency. One can employ different 
data fidelity and regularization terms such as those based on total variation (TV) and 
norms [45, 46]. Different selection can be effective for different problems. We refer interested 
reader to [45, 46] for detail numerical analysis and comparison of performances of different 
approaches. The primary focus here is the properties of the inverse transport problems in 
QPAT, not the details of the numerical implementation. 
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It is important to notice that the norm of the operator is in general smaU compared 
to that of and the identity operator X, especiaUy when the transport process is diffusive 
and the scattering is very isotropic. In those cases, [/(x, v) is roughly independent of v so 
that ||i^(t/)||Li(x) <^ 1- Thus {K{U)G)v is very small. In this case, when noise is present in 
the measured data, the reconstruction of the scattering coefficient is easily corrupted by 
noise. That was observed in the past in optical tomography [9, 10, 77]. We thus separate 
the reconstruction into the following two cases. 

3.2.1 Reconstructing (T, cr^) 

To reconstruct and T assuming is known, we solve the simplified least-squares problem 



mm 



/ I CI \ 



T 

Tan 



\1+P\\vl^^.^y\l. (3.17) 



We recover from the solution the coefficients T and a^. The solution of the least-squares 
problem is simply 



T 



1 f T 



X 



... C 



No 



V H ) 



(3.18) 

where the superscript * is used to denote the adjoint of an operator, and X* = X. The inver- 
sion of the regularized normal operator is achieved by a linear solver, not direct inversion. 

3.2.2 Reconstructing (a^, cr^) or (T, a^) 

To reconstruct either (a^, or (T, a^), we need to deal with an unbalanced system of 
equations caused by the smallness of the scattering component . We consider here only 
the case of reconstructing (a^, cr^). To obtain a similar algorithm to reconstruct (T, a^), 
we only need to replace the operator CJ^ in the algorithm by X. We solve the simplified 
least-squares problem 



mm 



To-., 



/ ^1 \ 



IP 



TcJa 



|2 
U2- 



(3.19) 



The normal equation for this minimization problem is 



i=l 



i=l 




E^g ra^ u 

i=l 'H 

E^g rs^-L 

i=l 'H 



(3.20) 



where Va and Vg are such that V = didig{Va^Vs). Instead of solving this normal equation 
directly to get a solution similar to (3.18), we perform Gauss elimination of the system to 
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obtain 



V C'-^'^ 



where 



i=l 1=1 i=l j=l 



and 



Ng Ng Ng Ng Ng 

K-^" = J2 ^*h, - [ A** A" ( E + p^a-^o) E A''* a] E A^^/^i- 

i=l 2=1 1=1 i=l i=l 



This motivates the following two-step reconstruction procedure. 



Step I. We reconstruct the coefficient Ta^ by solving the reduced linear system 

C^-'^Tas = /i""^. (3.22) 

In the construction of the operator C^^^ we need to apply the inverse of the operator 
^fji C^^^i + pT^a^a- We construct this inverse as accurate as possible. For large problems, 
we avoid constructing the inverse operator directly. Instead, we construct it implicitly. More 
precisely, to apply the inverse operator to any object z to get z = {J2i^=i ~^ P^aVa)~^ ^ ^ 

we solve the linear system ( ^fji + P^l^a)^ — z io maximum accuracy. 



Step II. Once the coefficient Tdg is reconstructed, we can eliminate it from the system 
to obtain a hnear system for the reconstruction of Tda- {^f=i^T^i + P^lPa)^^a = 

E5i>cr/^.-(E5iA^*A^)Ta,. 



3.3 Minimization-based nonlinear reconstruction scheme 

To solve the full nonlinear inverse problem, we reformulate the problem into a regularized 
least-squares formulation. More precisely, we minimize the following objective functional: 

1 

$(T,a„,0 = - J]||Ta„(n,)v-i/;||i2(a) (3.23) 

where Ui is the solution of the transport equation with the ith illumination, i.e., Ui solves 

V • VMj(x,v) + cr„(x)'Ui(x, v) = (Js(x)K(Mj) (x, v) inX .^24) 
'Ui(x,v) = ^i(x,v) on r_. 

The interior data collected for illumination gi is denoted by H*. 
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It is known that the functional (3.23) is Frechet differentiable with respect to the un- 
knowns provided that the coefficients satisfy the assumptions in (Al); see for instance [37, 
38, 79]. We can thus use gradient-based minimization technique to solve this problem. To 
calculate the Frechet derivatives of the objective functional, we use the method of adjoint 
equations [ ]. Let us denote by Wi the solution of the following adjoint transport equation 

-V •Vwi{-x,v) + aa{-x.)wi{-x,v) = as{-x.)K{wi){-x,v) + Taa{Taa{ui)^ - H^) in X 

i(;^(x, v) = on r+. 

(3.25) 

It is then straightforward to follow the standard calculations in [37, 38, 79] to show that the 
Frechet derivatives can be computed as follows. 

Theorem 3.5. The Frechet derivatives of ^ are given by 

^ar'^^ = 5]^(^^«^^^)v-^*)^«(«^)v>T)^.(^), (3.26) 

i=l 

,cra) = 2_^{{Taa{ui)^ - H*)T{ui)^ - {uiWi)^,aa)L^x), (3.27) 



1=1 



^s) = 2_^{{K{ui)wi)^, o's)l^{x) (3.28) 
i=i 

where (•, •)l2(x) denotes the usual inner product in L^(X). 

Note that due to the fact that the problem of reconstructing T is a linear inverse problem, 
the Frechet derivative of the objective functional with respect to T is independent of the 
adjoint solutions {wi\^^^. If we know and cfg but are interested in reconstructing T, then 
we can simply solve for T such that (1^? T) = for any test function T. This in turn gives 
us Taa{ui)v — i7* = 0, I < i < Ng^ which holds for every point x in Q. The least-squares 

solution of this overdetermined system is T = '}2i=i I Yl,i=i is^a{'^i)V}' Thus we 

need only solve the Ng transport equations and evaluate {ui)^ to reconstruct T. 

We use the limited memory version of the BFGS quasi-Newton method to solve the 
minimization problem. The details of the implementation are documented in [79, 77]. In 
our numerical experiments, we only attempt to reconstruct two of the three coefficients. We 
observe that the algorithm converges very fast, even from initial guesses that are relatively 
far from the true coefficients. This confirms the theory developed in the diffusive regime [1 4, 
16, 18] that is the problem is very well-conditioned when only two coefficients are sought. In 
the numerical simulation, we can add a small amount of regularization to the problem when 
the noise in the data is significant. This is done by adding a term of Tikhonov functional, 
such as ^||(T, Go) — (T, a^)!! (^1(^)^2 in the reconstruction of (T, a^) where T, a^, and are 
values obtained with a priori knowledge, in the objective functional (3.23). 
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4 The multi-spectral QPAT 



In practical applications of PAT, light of different wavelengths can be used to probe the 
properties of the medium at those wavelengths [25, 26, 27, 30, 33, 57, 66, 74, 76, 83, 84, 94, 
103, 105]. This is the idea of multi-spectral QPAT. 

Let us denote by A the set of wavelengths that can be accessed in practice, then the 
radiative transport equation in this case takes the form 

V • V2i(x, V, A) + cra(x, A)t^(x, V, A) = ^^(x, A)K(2i)(x, V, A) in X x A (il) 

2i(x, v,A) = 5'(x, v,A) on r_ X A, 

where the coefficients cra(x. A) and crs(x. A) are now functions of the wavelength. The interior 
datum collected in this case is now also a function of the wavelength 

//(x. A) = T(x, A) / a,(x, X)u{^, v, A)rfv. (4.2) 

Due to the fact that the equations for different wavelengths are de-coupled, we could not 
expect to reconstruct unknowns at one wavelength from data collected at other wavelengths, 
unless other a priori information is supplied. 

If the medium is non-scattering, the results in Section 2 enable us to reconstruct T(x, A) 
and cra(x. A) with two well-chosen illuminations 5'i(x, v. A) and 5^2 (x, v. A). The illuminations 
can be either collimated sources of the form g(x, A)(^(v — v') or point source in space of the 
form 0(v, A)(^(x — x'). We can simply perform wavelength by wavelength reconstructions 
using the methods developed in Section 2. 



4.1 Uniqueness of reconstruction 

The theory developed in [ ] for the diffusion equation confirms that we can reconstruct 
all three coefficients T, a^, and simultaneously with multi-spectral interior data under 
practically reasonable assumptions on the coefficients. Following [ ], we take the following 
standard model for the unknown coefficients: 

K 

a„(x, X) = J2 «fc(^)^a(x), <Ts{^, A) = /3(A)a,(x), T(x, A) = 7(A)T(x), (4.3) 

k=l 

where the functions {ak{X)}^=i^ /3(A) and 7(A) are assumed to be known. In other words, 
all three coefficient functions can be written as products of functions of different variables. 
Moreover, the absorption coefficient contains multiple components. This is the parame- 
ter model proposed in [30, 33, 57, 66, 94] to reconstruct chromophore concentrations from 
photoacoustic measurements. 

We have the following uniqueness result with the multi-spectral data. 

Corollary 4.1. Let (T,^^,^^) and (T,^^,^^) be two sets of coefficients of the form (4.3)^ 
satisfying the regularity assumptions in (Al). Assume that we have data from M(> K) 
different wavelengths {A^}^^;^ ^^^^ ^^^^ matrix cx {(^k{^m))^ l<k<K^l<m<M 
has rank K . Let A\ he the albedo operator that depends on wavelength. Then A\ = Ax 
tmpUes {T(x),KHx)}f^„a,(x)} = {t(x), KHx)}f^„ a,(x)}. 
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Proof. Corollary 3.2 implies that we can reconstruct uniquely 



^(x, A) = 7(A)T(x)a,(x, A), a(x, A) = a„(x, A) + /3(A)a,(x) (4.4) 
Take two wavelengths Ai and A2. We then have 

/^(x, Ai) ^ j{Xi)aa{-x, Xi) 

/i(x, A2) 7(A2)cra(x, A2) ' 

/5(A2)(t(x, Ai) - /3(Ai)a(x, A2) = /3(A2)a„(x, Ai) - /3(Ai)(t„(x, A2). (4.6) 

If /3(Ai)7(Ai)/i(x, A2) — /5(A2)7(A2)/i(x, Ai) 7^ at x G f^, we can solve the above system of 
equations to reconstruct aa(x, Ai) and aa(x, A2). We then obtain T(x), and ^^(x). Once 
we know T(x) and crs(x), we can reconstruct uniquely cra(x, A) for any A G A. Now select 
{Xm}m=i from A such that the matrix cx (Q^/c(Am)), I < k < K^l < m < M has rank K. 
We can reconstruct {o'^{'x.)}^^^ by solving the system X^^i <^fc(Am)<^a(^) = ^(x, A^), 1 < 
m < M. This completes the proof. □ 



4.2 Reconstruction methods 

The linearized and nonlinear reconstruction methods proposed in the previous section can 
be adapted to use the multi-spectral interior data. We would not repeat the whole algorithm 
again but just highlight the main modifications here for the linearized reconstruction with 
Born approximation. 

We can build an analogue of (3.11) 

(x) = J(T) + £«(Ta„) + CliTa,), (4.7) 



where I is the identity operator and the operators and are defined, respectively, as 

^li^^a) = ^+ /T(y)a<,(y,A)^^^(y;x)dy, (4.8) 

Cara,) = -a(A)^T(y)a,(y)i^^^^l^(y;x)dy. (4.9) 

Let us assume again that we collect the data for A^^ different illumination patterns and 
for each illumination patten, we have data for M different wavelengths {A^}^^^. We denote 
by Hi^x^^ 1 < i < Ng^l < m < M the ith data set at wavelength A^. The system of linear 
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integral equations with these Ng x M data sets can be written as 
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(4.10) 



Na^X 



J 



with X and C^- ^ 



being the evaluation of the operators defined in (4.8) and (4.9) respec- 
tively at source gi and wavelength A^. If we introduce the notation 



(Jn 



hi 



(a-«(x, Ai), . . . , (5-«(x, A^), . . . , a-a(x, Am))^, 

(^i,Ai5 • • • 5 ^i.Xm') • ' ' -) ^i.Xm) •) 



then this system is in the exact same form as (3.14). We solve this linear system in least- 
squares sense again in exactly the same ways as those presented in Section 3.2. Once 
(T(x), {aa(x, A^)}^^;l5 ^s(^)) reconstructed, we reconstruct the coefficient components 
{c^^}k=i by solving the linear system (in least-squares sense) locally (i.e. at each spatial 
location x G f^) 
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( ^^(x) \ 




a„(x, Ai) \ 
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1 (x) ; 
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a„(x, Am) / 



5 Numerical experiments 

We now present some numerical simulations with synthetic data to demonstrate the perfor- 
mance of the algorithms that we have presented in the previous sections. To simplify numer- 
ical computation, we consider only two-dimensional simulations, although the algorithms we 
have described are independent of spatial dimension. We will use both x and (x, y) to denote 
a point on the plane. The computational domain is the square ^ = (0, 2) x (0, 2). We denote 
hj ^} = dQ. We also denote by dQ\R, O^Ib, dfl\T the left, right, bottom and top 

sides of the boundary dQ respectively. For instance, dQ\L = {{x^y) | x = 0, G (0, 2)}. The 
unit sphere §^ is parameterized by an angle ^ G [0, 27r) so that the direction vector v can 
be represented as (cos^, sin0). As before, we will denote by {u)^ and {u)o the average of u 
on §1. 
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To generate synthetic data, we solve the transport problem on a spatial mesh that is four 
times as fine as the mesh used to solve the inverse problem. We then compute H in (1.3) 
by averaging in the four neighbor cells. This way, the interior data H contain automatically 
noise due to the mismatch in spatial discretization. However, we will still call the data 
constructed this way the "noiseless" data. In the case when the parameters that we are 
interested in are discontinuous in space, there is no way to recover the discontinuity in the 
coefficients exactly even with noiseless data. 

To get noisy data, we add random noise in the following way. Let H G (A^^^ being 
the total number of grid points) be a vector of noiseless data on the grid, then the noisy 
data vector H G is given by 

H = (J + eAA)H, AA = diag(Xi,...,X^J, (5.1) 

where Xj, j = 1, . . . , A^^^ are independent identically distributed Gaussian random variables 
with zero mean and unit variance, and e is the parameter that controls the noise level in the 
noisy data. For a particular value of the parameter e we say that the noise level is e • 100%. 
For example, for e = 0.05 we say the noise level is 5%. 



5.1 Reconstructions in non-scattering media 

In this section we present some numerical simulations in non-scattering media following 
the reconstruction methods presented in Section 2. To do the reconstruction, we cover the 
domain with 100 x 100 cells of uniform size whose nodes are given as 

= {xij = {xi,yj)\ Xi = iAx, Vj = jAy, ij = 0,1,..., 100}, 

with Ax = = 0.02. This is four times as fine as the grid used in next section for 
reconstructions in scattering media. 




5.1.1 Reconstructing 

We show now some reconstructions of the absorption coefficient. 




Figure 3: Reconstructions of a piecewise constant absorption coefficient with a collimated 
source. Left to right: true absorption coefficient a^, interior data H, a a reconstructed with 
noiseless data and reconstructed with noisy data (noise level 5%). 

In the first numerical simulation in this group, we perform a reconstruction of a piecewise 
constant absorption function using a collimated source located on the left side of the bound- 
ary pointing inside the domain. The source is ^'(x, v) = Xdn\L^{'^ ~ ^0 with = (1^0). 
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The absorption function consists of a background a a = 0.1 cm ^ and three disk inclu- 
sions f^i = {x G I |x - (0.6,0.6)1 < 0.2}, = {x G I |x - (1.4,0.6)| < 0.2} and 
= {x G I |x — (1.0,1.5)1 < 0.3} with values (Ja\Q^ = 0.3 cm~^, (Ja\Q^ = 0.3 cm~-^ 
and (Ja\Q^ = 0.2 cm~^ respectively. The reconstruction results with noiseless and noisy data 
are presented in Fig. 3. The reconstruction is almost perfect when noiseless data is used. 
When noisy data (e = 0.05) is used, we can clearly see a degeneration of the quality of the 
reconstruction. However, the degeneration is very small, comparable to the noise level of 
the data. 
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Figure 4: Reconstructions of a piecewise constant absorption coefficient with cone-limited 
source. Left to right: true absorption coefficient a^, interior data i7, a a reconstructed with 
noiseless data and reconstructed with noisy data (noise level 5%). 

Next we perform two reconstructions using a cone-limited source. The setup is as de- 
picted in Fig. 2. The boundary condition g{'x^9) corresponds to a uniform isotropic line 
source of unit intensity concentrated on a segment {(x,^) | x G (0,2),^ = —2} which re- 
sults in 5'(x, ^) being non-zero only on dQ\B- The half-aperture angle for such source is 
7r/4, so we only keep track of the solution u{x^y^6) for 6 G [7r/4, 37r/4]. This segment is 
uniformly discretized with 50 nodes to generate the data and with 40 nodes to solve the 
inverse problem using the algorithm from Section 2.3. 
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Figure 5: Reconstructions of a smooth absorption coefficient with cone-limited source. Left 
to right: true absorption coefficient a^, interior data i7, a a reconstructed with noiseless data 
and da reconstructed with noisy data (noise level 5%). 
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We show in Fig. 4 the reconstruction of a piecewise-constant absorption coefficient with 
the cone-hmited source. The absorption coefficient (in units of cm~^) is given by 

a„(x) = 0.1 + 0.1xQi(x) + 0.2xQ2(x) + 0.3xn3(x), (5.2) 

with the rectangular inclusion Qi = [0.4,1.2] x [0.2,0.8], the smaller disk inclusion Q2 = 
{x G I |x - (1.6, 1.0)1 < 0.3} and the larger disk inclusion = {x G | |x - (0.6, 1.4) | < 
0.4}. Thus the absorption coefficient in the background is = 0.1 cm~^ while that in 
the three inclusions it takes the values aa\n^ = 0.2 cm~-^, aa\n2 = 0-3 cm~-^ and cFaln^ = 
0.4 cm~^ respectively. The quality of the reconstruction is comparable to that in the previous 
numerical experiment in Fig. 3. Smoother absorption coefficients can be reconstructed with 
similar quality. To demonstrate that, we show in Fig. 5 the reconstruction of the absorption 
coefficient that is a sum of a Gaussian and linear functions given by 

aa{x, y) = Ax + By + C + De {i-^-^fHy-^-^f) / si ^5^3) 

where the parameters ^4, 5, C, and Sd are chosen so that 0.1 cm~^ ^ era ^ 0.3 cm~^. We 
performed reconstructions for many different choices of smooth absorption coefficients. The 
quality of the reconstruction in Fig. 5 is representative of the typical reconstruction quality 
that we obtained in those experiments. 




Figure 6: Relative error 8(j in the reconstruction of the absorption coefficient with a cone- 
limited source versus noise level £h in the data. Noise levels from to 30%. Left: recon- 
struction of piecewise constant absorption coefficient (5.2); Right: reconstruction of smooth 
absorption coefficient (5.3). Perfect linear stability £cr = £h is shown as a dashed line for 
reference. 

To characterize the stability of the reconstruction more precisely, we plot in Fig. 6 the 
relative error S^- in the reconstruction of the absorption coefficients (5.2) and (5.3) versus 
the noise level £h in the data used, with S^r and £h defined respectively as 

\\d-a - CTa\\l2 ||H-H||/2 

= M M ' = — liTTTl (5-4) 
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where a a G and a a ^ 1^^^ are the true and reconstructed absorption vectors respec- 
tively. NumericaUy the method appears to have hnear stabihty, with the piecewise constant 
case being shghtly worse than the smooth one. This is typicaUy due to an imperfect resolu- 
tion of the boundaries of the inclusions. Also, there is some residual error in the noiseless 
(e = 0) case due to the mismatch between the fine grid and the reconstruction grid. 

5.1.2 Reconstructing (T, a^) 

To reconstruct both the Griineisen coefficient and the absorption coefficient, we use data 
collected from two coUimated sources located on the left and right sides of the boundary 
respectively, ^i(x, v) = XdQ.\i^K^ - v') and ^2(x, v) = Xa^7|i^^(v + v') with = (1, 0). 




Figure 7: Reconstructions of (top row) and T (bottom row) using a pair of coUimated 
sources. Left to right: true coefficients, coefficients reconstructed with noiseless data and 
coefficients reconstructed with noisy data. The noisy data contains 5% random noise (e = 
0.05). 

The background absorption coefficient is = 0.1 cm~^ and the background Griineisen 
coeflBcient is T = 0.5. There are four inclusions, two for the absorption coefficient located in 
f^i = {x G I |x- (0.5,0.5)1 < 0.3} and = {x G | |x - (1.5, 1. 5) | < 0.2} respectively 
and two for the Griineisen coefficient located in = {x G | |x — (0.5, 1.5) | < 0.2} and 

= {x G I |x — (1.5,0.5)1 < 0.3} respectively. The coefficients inside the inclusions 
are Ga\vL-Y = 0.2 cm~-^, Ga\vL2 = 0.3 cm~-^, T|^3 = 0.6 and Tl^^ = 0.7 respectively. We 
perform the reconstruction with both noiseless data and noisy data polluted with 5% additive 
random noise. The reconstruction results are presented in Fig. 7. Other than the phantom 
inclusions in the reconstructed Griineisen coefficient at the locations of the inclusions of the 
absorption coefficient, the quality of the reconstructions is very high, comparable to the 
previous reconstructions in the cases of one unknown coefficient. The phantom inclusions in 
the reconstructed Griineisen coefficients are caused by the inaccuracy of the reconstruction 
of the absorption coefficient at the boundary of square inclusions. This inaccuracy originates 
from the differentiation of the quantity In ^ which contains noise coming from mismatch 
between the forward and inversion grids. For the reconstruction from noisy data, the noise 
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added to the synthetic data in random but have only low frequency components. If the data 
contain very high frequency components as in the previous cases, numerical differentiation 
of the quantity In ^ in the algorithm would yield even larger noise in the reconstruction 
that would bury the true coefficients. Averaging from multiple reconstructions would be 
necessary to get a clean image. We would not address this issue in detail now but leave it 
to future study. 

5.2 Reconstructions in scattering media 

Here we present some numerical simulations for scattering media, i.e. the case when 7^ 0. 
In this case, we do not have analytical reconstruction formulas to work with. We thus use 
the linearized reconstruction and minimization-based nonlinear reconstruction schemes. 

5.2.1 Numerical setup 

The transport equations are solved with a finite volume scheme for the spatial variable and 
a discrete ordinate method for the angular variable. The domain is covered by 50 x 50 cells 
of uniform size whose nodes are given by 

Qh = {^ij = {xi,yj)\ Xi = iAx, yj = jAy, ij = 0, 1, . . . ,50}, 

with Ax = Ay = 0.04. We discretize §^ into 128 uniformly distributed directions with 
identical quadrature weight: 

^io = {v, : V, = (fc - 1) * A^, = 1, . . . , 128}, 

where = 27r/128. Note that this discretization is not as fine as that used in the previous 
case when scattering is not presented in the problem. This is only done here to save com- 
putational cost. It should not be regarded as a limitation of the reconstruction algorithms 
for scattering media. 

The scattering kernel is chosen as the Henyey-Greenstein phase function [9, 54, 97] 

where 77 G [0, 1] is the anisotropy factor, which measures how peaked forward the phase 
function is. The larger rj is, the more forward the scattering. The anisotropy factor is often 
used to define the so-called effective scattering coefficient through = (1 — 77)0-5. 

For more details on the forward solver, we refer to our previous publications [78, 79]. 

5.2.2 Reconstructing 

We first reconstruct the absorption coefficient, assuming that both the scattering coefficient 
and the Griineisen coefficients are known. The setup is as follows. The medium consists 
of a homogeneous background absorbing medium with absorption coefficient = 0.2 cm~^ 
and two absorbing inclusions. The first inclusion occupies = [0.4, 0.8] x [0.4, 0.8] with the 
absorption coefficient = 0.3 cm~^. The second inclusion occupies = [1-2, 1.6] x [1.2, 1.6] 
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Figure 8: Reconstructions of the absorption coefficient. Top row: the data H = Taa{u)^ 
(left), aa reconstructed using Born approximation (middle) and reconstructed using non- 
linear iteration (right) in an isotropic medium; Bottom row: same as the top row but in an 
anisotropic medium. 

with the absorption coefficient = 0.1 cm~^. We first performed two reconstructions with 
the linearization method: a reconstruction in an isotropically scattering medium with rj = 0^ 
as = 8 cm~^ and a reconstruction in an anisotropically scattering medium with rj = 0.9, 
as = 80 cm~^. The Griineisen coefficient is always T = 0.5. The synthetic data used in 
these reconstructions is noiseless in the sense described at the beginning of this section. The 
results of the reconstructions are shown in Fig. 8. The quality of the reconstructions is very 
high (although shown on a coarser grid than those in the previous section). Note that the 
fast decay of field from the line source makes the second inclusion hardly visible in the data 
H = Taa{u)v plots. This is one of the main reason that the quantitative step of PAT is 
needed. The reconstructions using the nonlinear reconstruction algorithm are presented in 
the right column of Fig. 8. 

5.2.3 Reconstructing (T, a^) 

We now perform numerical simulations where we reconstruct both the Griineisen coefficient 
and the absorption coefficient. The scattering coefficient is fixed at = 8 cm~^ and the 
anisotropic factor ij = 0. The background absorption coefficient is = 0.1 cm~^ and 
the background Griineisen coefficient is T = 0.5. There are four inclusions, two for the 
absorption coefficient located in = [0.4,1.2] x [0.4,0.8] and ^2 = [0.4,0.8] x [1.2,1.6] 
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Figure 9: Reconstructions of the Griineisen and the absorption coefficients. Top left: true 
coefficients (T, cr^) Top right: two data sets (iifi, H2) used in the reconstruction; Bottom 
left: reconstructed (T, a^) with two data sets; Bottom right: reconstructed (T, a^) with 
eight data sets. 

respectively and two for the Griineisen coefficient located in = [0.4, 1.6] x [1.2, 1.6] and 
Q4 = [1.2,1.6] X [0.6,1.0] respectively. The coefficients inside the inclusions are aa\n^ = 
0.2 cm~-^, aa\n2 = 0.3 cm~-^, T|^3 = 0.6 and Tl^^ = 0.7 respectively. We perform the 
reconstruction with data polluted with 5% additive random noise added the same way as 
in the previous section. The coefficients (T^aa) reconstructed with two data sets generated 
with sources 5'i(x,v) = XoQl 5^2 (x,v) = XdQn are shown in Fig. 9 (the left two plots 
in the bottom row). The quality of the reconstruction is quite high again even though it is 
slightly lower than that in the reconstructions in the non-scattering regime, such as those 
shown in Fig. 7. The quality of the reconstruction can be improved by averaging out noise in 
the data through the usage of additional data sets. This is demonstrated in Fig. 9 (the right 
two plots in the bottom row) where we show the reconstruction of the coefficients (T,^^) 
using eight data sets. The eight sources used are the rotation of the source 5'i(x) = XoQl 
and g2 = XdQL^i'^ ~ (I5 0)) ^^^^ ^/^^ ^r, and 37r/2 around the center of the square domain 

5.2.4 Reconstructing (a^, (Js) 

We now fix the Griineisen coefficient T = 1 and attempt to reconstruct both the absorption 
and the scattering coefficient. We show in Fig. 10 the reconstructions of the objective 
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Figure 10: Reconstructed absorption and scattering coefficients given in (5.6). Top row: (a^, 
as) reconstructed with noiseless (left two plots) and noisy (right two plots) data; Bottom 
row: the corresponding relative differences between reconstructed and real coefficients. 



absorption and scattering coefficients: 

(Ja{x,y) = 0.2 + 0.1sin(7rx), crs{x,y) = 8.0 + 2.0sin(7r?/). (5.6) 

Shown are reconstructed {(Takers) (top row) with eight noiseless and noisy data sets and the 
corresponding pointwise relative error (bottom row), defined for (resp. a^) as ^^^^^ 
(resp. ^^^), in the reconstructions. Except for the relatively large error on some part of 
the boundary, the error in the reconstructions is comparable to the noise level in the data. 
This is also confirmed in the reconstruction of piecewise constant coefficients such as those 
shown in Fig. 11. The piecewise constant absorption coefficient but the smooth scattering 
coefficient are given as 

(Ja{x,y) = 0.2 + 0.1sin(7rx), crs{x,y) = 8.0 + 2.0 sin(7r2/) (5.7) 

where the first absorbing inclusion is located in Qi = [0.4,0.8] x [0.8, 1.2] and the second 
absorbing inclusion is located in Q2 = [1-4, 1.6] x [0.2, 1.8]. For piecewise constant coefficient, 
the reconstruction error occurs on the boundary of the inclusions. This is true also in the 
noiseless data case due to the fact that the synthetic data is generated by averaging quantities 
on a finer mesh. The Tikhonov regularization we employed in the numerical schemes also 
contribute to the smoothing on the boundary of the inclusions. 
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Figure 11: Reconstructed absorption and scattering coefficients given in (5.7). Top row: (a^, 
as) reconstructed with noiseless (left two plots) and noisy (right two plots) data; Bottom 
row: the corresponding relative differences between reconstructed and real coefficients. 

5.2.5 Reconstructing (T, a^, cTs) with mult i- spectral data 

The last numerical simulation is devoted to the simultaneous reconstructions of all three 
coefficients with multi-spectral interior data. The coefficients take the forms given in (4.3). 
To simplify the presentation, we consider only an absorption coefficient that has two com- 
ponents, i.e. K = 2. The spectral components of the coefficients are given as follows: 



ai{X) = ^, a2{X) 



3/2 



7(A) = 1, 



where the normalization wavelength Aq is included to control the amplitude of coefficients. 
These weight functions are by no means what they should exactly be in practical applica- 
tions. However, the specific forms do not have impact on the results of the reconstruction. 
The spatial components of the coefficients are given as 



T(x) = 0.5 + 0.4tanh(4x - 4), 
ai(x) =0.1 + 0.1x^7, +0.2x^^,, 



as{'x.) = 8.0 + 2.0 sin(7rx) s'm^ny) 
a2(x) = 0.2 -0.1x^^3 +0.1x^74 



(5.9) 



where Qi = [0.4,1.6] x [0.4,0.6], = [0.4,0.8] x [0.8,1.2], = [0.4,1.6] x [1.4,1.6], and 
= [1.2,1.6] X [0.8,1.2]. The anisotropic factor r] = 0. We performanced numerical 
reconstructions using four wavelength-dependent sources. For each source, we have data 
for four different wavelengths. The results of the reconstructions using noiseless data are 
presented in Fig. 12 and those using noisy data are shown in Fig. 13. We observe similar 
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Figure 12: Simultaneous reconstruction of the absorption and the scattering and the 
Griineisen coefficients with multi-spectral data. Top row: reconstructed T, a^, al and 
a^; Bottom row: the corresponding pointwise relative error in the reconstructions. The data 
used in this simulation are noiseless. 



reconstruction qualities (as can be seen in the plots of the pointwise relative error) to the 
two-coefficient cases in the previous sections. 

6 Concluding remarks 

We studied the quantitative photoacoustic tomography problem with the radiative transport 
model, aiming at reconstructing multiple physical coefficients simultaneously using the data 
collected from multiple illuminations. We showed that in non-scattering absorbing media, 
we can reconstruct both the absorption and the Griineisen coefficients simultaneously in a 
stable manner using only two sets of interior data. Moreover, in this case we derived explicit 
reconstruction formulas for the problem with particular choices of illuminations (coUimated, 
point and cone-limited sources). In scattering media, we show, based on the result in [^^]^ 
that one can reconstruct two of the absorption, the scattering and the Griineisen coefficients 
stably when more data, i.e., data given by the full albedo operator, is available. To re- 
construct all three coefficients simultaneously, we proposed to use interior data collected at 
different optical wavelengths. We show that with some realistic a priori knowledge, mainly 
the knowledge of the spectral denpence of the coefficients, we can reconstruct stably the 
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Figure 13: Same as in Fig. 12 except that the data used contain 5% random noise. 

spatial component of aU three coefficients simultaneously. 

Besides the analytical reconstruction strategy for the non-scattering problem, we pro- 
posed a linearized reconstruction method based on Born approximation to the original in- 
verse problem as well as a nonlinear reconstruction method based on numerical minimization 
techniques for QPAT for scattering media. We show numerically that the reconstruction is 
very stable. In fact, our numerical experiments showed that, assuming that the data mea- 
sured with ultrasound is accurate enough, the nonlinear least-square formulation of the 
inverse transport problems with interior data problem behaves like a convex optimization 
problem and thus can be solved efficiently and accurately; see the numerical reconstructions 
in Section 5. 

There are many important issues in quantitative PAT that need to be addressed in the 
future. For instance, in the theory developed in [14, 18, li ] for the diffusion model, one can 
reconstruct both the absorption and the scattering coefficients with only two "well-chosen" 
illuminations, with little a priori assumptions on the coefficients. It would be interesting to 
see how to generalize that to the radiative transport model (1.1) with 7^ 0. An equally 
important question is whether or not we can derive any analytical reconstruction procedure, 
similar to those developed in diffusion-based theory [18, 14, 16], for the inverse transport 
problem in scattering media. 

On the numerical side, when the noise present in the data is significant, we need to 
regularize the reconstruction. The regularization we adopt here is the usual Tikhonov reg- 
ularization which yields smooth solutions among all possibilities. In certain applications. 
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we might know a priori that the coefficients to be reconstructed are piecewise constant, 
such as those presented in some of the numerical simulations in Section 5. In this case, 
alternative regularization strategies such as the total variation (TV) regularization might be 
more appropriate. In [14, 47^ 4^]^ regularization for non-smooth coefficients have been 
considered in the diffusion case. It would be interesting to see the performance of that in 
the transport case. We plan to explore this issue in the future. 
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